Role of the self-interaction error in studying first principles chemisorption on graphene 
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Adsorption of gaseous species, and in particular of hydrogen atoms, on graphene is an important 
process for the chemistry of this material. At the equilibrium geometry, the H atom is covalently 
bonded to a carbon that puckers out from the surface plane. Nevertheless the flat graphene geometry 
becomes important when considering the full sticking dynamics. Here we show how GGA-DFT 
predicts a wrong spin state for this geometry, namely S z =0 for a single H atom on graphene. We 
show how this is caused by the self-interaction error since the system shows fractional electron 
occupations in the two bands closest to the Fermi energy. It is demonstrated how the use of hybrid 
functionals or the GGA+?7 method an be used to retrieve the correct spin solution although the 
latter gives an incorrect potential energy curve. 



I. INTRODUCTION 

Thanks to its peculiar semi-metallic band structure, 
graphene is a very promising material for future silicon- 
free nano electronics [3, 0- In particular, graphene's 
charge-carriers mobility is indeed extraordinarily high, 
with mean- free paths in the order of micronsjl, Qj. How- 
ever, for the fabrication of logic devices the absence of a 
band gap is a major issue since it does not allow a com- 
plete turn-off of the current, hence high on-off ratios [HQ. 
One straightforward possibility for graphene band-gap 
engineering is to use the adsorption of radicals or small 
molecules to create 7r defects, breaking the equivalence 
of its two sublattices. In this perspective several studies 
have been performed on isolated^ 0|, clusters and 
even superlattices of adsorbates bonded on graphene, and 
most of them are based on the density functional theory 
(DFT). Among these hydrogen atoms adsorption is by 
far the most studied case, also because of the many im- 
plications in other fields (l2l. [l3|. 

Hydrogen chemisorption dynamics on graphene is not yet 
completely understood. In the adiabatic picture when a 
H atom impinges on graphene and binds at its top site 
it induces one carbon atom to move out from the plane, 
"puckering" the surface. On the other hand, if the incom- 
ing species moves fast enough toward the graphene layer, 
sticking can occur faster than surface reconstruction; in 
this case the substrate can be considered as rigid. The 
planar geometry may thus play an importance role in 
the adsorption dynamics, and the accuracy of DFT cal- 
culations in this case is crucial in order to build reliable 
potential energy surfaces for quantum dynamics studies. 
In this article we analyze the spin properties of the sub- 
strate when a H atom is forced to chemisorb on a flat 
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graphene sheet. We show how the loss of spin polariza- 
tion is a fictitious feature of semi-local generalized gra- 
dient approximation (GGA) functionals due to the self- 
interaction error (SIE). Then we will show how the cor- 
rect ground state magnetization can be achieved using 
hybrid functional or GGA+U approaches, even though 
the latter approach fails in reproducing the correct H- 
graphene potential energy curve. 

II. COMPUTATIONAL METHODS 

In brief, periodic density functional theory as imple- 
mented in the VASP package [lil Il5j has been used 
throughout. A GGA-PBE functional was used and the 
plane wave basis set was limited to a 500 eV energy cut- 
off. For the inner electrons we rely on the frozen core 
approximation using PAW pseudo-potentials [l6j . 
The reciprocal space was sampled by T centered k-point 
grids, whose meshes were chosen depending on the su- 
percell size, in any case never more sparse than 6x6x1. 
The graphene unit supercells used here range from a 2x2 
to a 5x5: all of them have a vacuum region along the c 
axis of 20 A in order to guarantee a vanishing interaction 
between periodically repeated images. 
It has been shown recently that for hydrogen adsorption a 
5x5 supercell is still not bigenough to extract adsorption 
energies at meV accuracy [8|, although this goes beyond 
the aim of this work. Further details about the compu- 
tational setup can be found in ref.[8j|. 

III. RESULTS AND DISCUSSION 

When a radical species chemisorbs onto graphene or 
graphite surfaces the most favorable outcome is the for- 
mation of a covalent bond with one of the surface car- 
bons, i.e. at a top site. The simplest case of radical 
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is a single (neutral) hydrogen atom. Applying Valence 
Bond (VB) arguments to this reaction one can predict 
that as the atom approaches the substrate plane it in- 
teracts with a 7r electron of graphene, triggering orbital 
re-hybridization of the C atom from a planar sp 2 to a par- 
tially tetrahedral sp 3 configuration^, 0|. At long range 
the adsorbate - substrate interaction is purely repulsive 
since there is no unpaired electron on graphene available 
to bind hydrogen. At short range however a low-lying 
spin-excited state in which two tt electrons lying on oppo- 
site, non-overlapping ends of a benzene ring would give 
rise to an attractive, barrierless interaction with the H 
Is lone electron. Hence, an avoided crossing between 
these two doublet curves occurs giving rise to an activa- 
tion barrier to chemisorption (see Fig. 5 in ref.Q). When 
the graphene sheet is allowed to reconstruct this process 
reads as a sp 2 — sp 3 orbital re- hybridization of the carbon 
involved in the bonding that turns partially tetrahedral, 
puckering out 0.6 A from the layer plane. 
The graphene lattice is a bipartite system, made of two 
equivalent sublattices, each made of every second carbon 
atom. The equivalence of the two honeycomb sublattices 
is responsible for the particle-hole symmetry in graphene, 
and for the peculiar conical intersection at Ep between 
the valence and conduction bands [13] ■ A chemisorbed 
species creates a defect in the aromatic network, hence 
an imbalance between the number of occupied sites of the 
two 7r sublattices (jia and ns respectively). According 
to a theorem formulated by Inui et. al. within tight- 
binding theory, whenever an imbalance (vacancy) is in- 
troduced in a bipartite lattice this gives rise to \ua — 
zero-energy midgap states, localized on one sublattice 
only flal - Moreover, following the second Lieb's theo- 
rem since the total number of electrons is odd the 
total magnetization for non-metallic systems has to be 
S 2 = \n A -n B \/2- Thus, for a single defect S 2 = l/2, or 1 
Ms- 

DFT calculations confirm this picture: the hydrogen 
atom introduces a defect in one of the two sublattices, 
breaking one among the many aromatic bonds around 
the tetrahedral carbon. This implies that an unpaired 
electron can be delocalized by a "bond switching" process 
along the other sublattice, made of every second carbon 
atom. In the energy spectrum this reads as a flat band 
at the Fermi level i.e. in a midgap state occupied by one 
single spin projection only[§|, [20 . 

According to previous studies [a 3 we found at the GGA- 
DFT level that keeping the substrate planar thwarts the 
sp 2 — sp 3 re-hybridization; this is enough to weaken the 
attractive C-H interaction, but also to corrupt the sys- 
tem's aromatic character. Nevertheless, because the VB 
arguments concerning the crossing of two spin states 
hold, only a meta-stable C-H bond can form. 
From our DFT calculations we notice how the total spin 
for H adsorbed on flat graphene is lower than expected 
at the local minimum geometry. In FigfT] is shown the 
value of magnetization (left panel) along the adsorption 
path for several surface coverages, together with the total 
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Figure 1: Left panel: total spin vs. distance from the surface 
for a H atom Adsorbing on flat graphene. The coverages are 
the following: 9=0.031 (full line), 0.055 (dashed) and 0.125 
ML (dotted). 

Right panel: adsorption potential curves for H (full line) to- 
gether with the same curves for the spin unpolarized case 
(squares) and for the fixed magnetization 1 (j,b (circles). The 
full line shows the adiabatic (non-constrained) curve. Zero 
energy here is set asymptotically away from graphene. 



energies for the spin-polarized and unpolarized solutions 
(right panel) . When the radical is far from the surface its 
magnetic moment is correctly 1 \ib'- this corresponds to 
an electron lying in the H Is orbital, while graphene elec- 
tronic structure remains intact. As the atom approaches 
the graphene layer along the normal direction, at a given 
critical height from the surface (z c (H) ~ 1.25 A) the sys- 
tem's spin drops. The minimum value for the total spin 
depends upon the coverage. When coverage is low enough 
the spin is eventually quenched down to zero. Pushing 
the adsorbate closer to the carbon atom the magnetiza- 
tion tends to increase again. We tested that the same 
picture also holds for other small organic radicals with 
slightly different critical heights, weakly dependent upon 
the supercell size (coverage). 
When comparing the adsorption curves computed with 
a magnetization fixed to and 1 [is as in FigfTJ (right 
panel), it is clear that the non-magnetic (non-polarized) 
solution becomes more stable than the magnetic one be- 
yond the critical point z c . 

To a closer look (see FigJ^b) we notice that the occu- 
pied hydrogen s orbital and its (empty) affinity level get 
closer in energy when approaching graphene consistently 
with the Newns- Anderson model |2ll . |22j . Then around 
z c they approach Ep, become degenerate and equally 
occupied by a fractional number of electrons. In the sys- 
tem's density of states, a gap opens at the point of the 
spin quenching together with the formation of a broad 
partially occupied peak at the Fermi energy, symmetric 
for both spin projections. 

If the host were a metal, then spin-flip scattering between 
conduction electrons and a magnetic impurity might lead 
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Figure 2: a: Occupancy for the spin up (circles) and spin 
down (squares) bands closest to Ef along the whole adsorp- 
tion path. Dashed line: total magnetization (in ^b)- b: 
Eigenvalues for the bands shown above, the dotted line rep- 
resents the value of the Fermi energy. 



to a singlet ground state (Kondo effect). However in our 
case there is not a delocalized free-electron-gas-like sur- 
face state (a Shockley state) that is free to screen the 
impurity. Moreover the s electron does not belong to a 
localized orbital decoupled from the substrate such as for 
d metals. The Kondo effect is a many body feature that is 
not well described by DFT. Anyway this technique can 
describe correctly the spin quenching on metallic sub- 
strates, such as for a hydrogen atom impinging on Cu, 
Ag or Al(lll) surfaces where this s pin transition is a sig- 
nature of non-adiabatic effects [23L l24j . 
We found that radical adsorption on planar graphene is 
a situation in which the (many-electron) self-interaction 
error (thereafter SIE, also known as delocalization 
error 12 51) is particularly severe. Following the notation 
in ref[26] this is a chemical reaction with two centers 
(m = 2) and one electron (n = 1), since none of the 
electrons of graphene may be directly involved in bond- 
ing if they cannot re-hybridize to a tetrahedral sp 3 state. 
Systems with a fractional n/m ratio are also known to 
exhibit of large static electronic correlation [27} ■ 
A well known system with n/m — 1/2 is the molec- 
ular ion dissociation. When the two atoms are far apart 
both the local density approximation (LDA) and GGA lo- 
cal functionals give as the most stable ground state half of 
an electron on each of the two degenerate atomic orbitals, 
which is physically not correct. Asymptotically the or- 
bitals on the two fragments are degenerate, so this frac- 
tional occupation solution should be degenerate with any 
other possible electronic arrangement such as one filled 
and one empty orbital [28[ . Fractional charge (and spin) 



on the two fragments is due to the over-delocalization 
that is a direct consequences of the SIE. Indeed the LDA 
and GGAs functionals are designed to correctly repro- 
duce the system's total density and the on-top pair den- 
sity, but fail in pathological systems to reproduce spin 
densities because of the SIE [5^, [3(j. This non- integer 
orbital occupation is directly related to the break down 
of sum rules over the exchange hole density n x (r, r') 

/ ^(r,r')dr'=^/ a , ff ^ (1) 
J ^ n(r) 

that when the orbital occupation / becomes fractional 
sums to a value within larger than the correct -1. Thus 
local functionals give an energy whose behaviour is con- 
vex for a fractional number of electrons instead of being 
linear according to Janak theorem [3l|. For this reason 
for open systems a delocalized situation with fractional 
charge turns out to be more favored [25]. 
For H adsorbed on flat graphene the spin quenching is 
similarly due to a fractional spin situation: the "split- 
ting" of one electron in two different degenerate bands 
(originating from H s and C p z ), with opposite spin pro- 
jections (see Fig|3]). Here the occupation number in these 
bands can fluctuate, a sign that SIE is particularly severe 
[32j. This picture is confirmed by the convex behaviour 
of the system energy for fractional band occupation as 
shown in FigOl obtained by constraining a given occu- 
pation within the two bands. Note that two degenerate 
orbitals every arrangement of electrons (even for frac- 
tional numbers) should be perfectly degenerate. In this 
case however, the energy minimum lays exactly at the 
unpolarized solution: 0.5 occupancy of the two bands. 

Another indication of the SIE is the worsening of the 
spin quenching at low coverages, hence for larger and 
larger super-cells, as shown before in FigJT] With the 
SIE being cause of the delocalization, the fractional oc- 
cupancy might not be maximal when the unit cell is not 
large enough to accommodate all the delocalized electron 
density [33j. 

A major difference with the H^~ prototype case is that 
here there is no fractional spin asymptotically for H since 
the involved orbitals here lay far below graphene Fermi 
energy. For other and more electronegative monovalent 
species such as F and OH fractional charges appear also 
asymptotically, similarly to the case of dissociations of 
heteronuclear diatomics [3^ |. 

To prove further that the failure in representing the to- 
tal spin of the system comes from the approximate nature 
of the density functionals we tested the performances of 
a (non-local) hybrid functional. Hybrid functionals com- 
bine the GGA exchange and correlation term (convex for 
fractional electron number) with the Hartree-Fock (HF) 
"exact" exchange term. Since HF energies have instead 
a concave behaviour for fractional electron numbers they 
often yield over-corrected results [33] . For this reason the 
use of a fraction (usually one fourth) of exact exchange to 
correct the DFT E xc functional gives much better results, 
at least for non- metallic systems [35| . We then chose the 
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Figure 3: Calculated total energy vs. magnetization for one 
H atom adsorbed on flat graphene. When the magnetization 
is -f or f fiB one of the two degenerate bands is occupied 
while the other is empty. For the non-polarized, /is, case 
both bands are occupied by half of an electron each. The full 
line represents the correct degenerate behaviour for fractional 
electron numbers; filled circles are the GGA-DFT results. 




Figure 4: Density of states for the H-flat graphene system ob- 
tained for a 2x2 supercell. PBE results for: a) clean graphene, 
b) H adsorbed on reconstructed (puckered) graphene, c) H 
adsorbed on flat graphene (M=0/is). Hybrid-PBEO results 
(1(j,b) are shown in d) for comparison. As a guide to the eye 
the Fermi energy is shown as dashed vertical line. 



PBE0 functional, which mixes 25% of HF with 75% of 
PBE exchange, and uses the full PBE correlation [36]: 



pPBEQ _ rpPBE 

xc xc 



0.25(£f F - E™*) 



PBE\ 



Hybrid functionals are orbital dependent, i.e. non-local 
in space: this is a major issue when employing plane 
wave codes where the number of orbitals depends upon 
the supercell volume. The computational effort needed 
for this kind of calculations is thus much larger com- 
pared to GGA. Hence we could not reach low coverages, 
restricting our system to a 2x2 supercell (0=0.125 ML). 
Within hybrid functional DFT the system's total spin for 
H on flat graphene layer is 1 [Ib (<Sz=1/2), and the asso- 
ciated spin density is correctly localized on one sublattice 
only. A comparison between the density of states com- 
puted with PBE and PBE0 is shown in FigOfl GGA can 
represent pretty well the adiabatic chemisorption mech- 
anism, namely the band gap opening and the zero en- 
ergy midgap state that splits by the exchange interac- 
tion into a filled and an empty states with opposite spin 
(Fi^jUa and b). For the flat geometry, GGA gives a broad 
feature straddling across the Fermi energy (Fig|4]c). For 
larger supercells (at lower coverages) the peak becomes 
fully symmetric for the two spin projections, giving rise 
to the un-polarized state. On the other hand PBE0 can 
reproduce well the occupied midgap state. As it is widely 
known standard DFT tend to underestimates the band 
gaps, again due to the SIE[37|: for this reason the PBE0 
band-gap is about 50% larger than that of PBE. 
It has also been proposed that an "on site" repulsion term 
such as in the LSDA+{7 [38] approach can help to control 
the SIE in case of fractional occupation [39| . The on site 
Coulomb term U acts as a penalty for the occupation 
of the two degenerate bands at the Fermi energy, and 
can thus reproduce the correct ground state. Note that 
LSDA+J7 and GGA(PBE)+C7, as implemented here, give 
practically identical results in this case. This approach 
is much less computationally expensive than the hybrid 
functionals, so we could study the full adsorption path. 
As for PBEO the total magnetization is 1 /is at every C-H 
distance, and the spin density correctly is localized either 
on the H atom or on one graphene sublattice. From our 
tests a Coulomb interaction of 15 eV was enough to re- 
trieve the correct ground state spin: a value not far from 
20.08 eV, already successfully used to describe carbon it 
electrons in similar approaches jiol l4l| . 

We would like to stress that the GGA+U approach is 
not a rigorous way to avoid self interaction, and hence 
some care has to be taken when interpreting these re- 
sults. While it is relatively easy to predict which has to 
be the correct total spin of the system upon physical ar- 
guments, it is more challenging to judge the quality of 
total energies. Indeed, it can be seen from Fig(S] that 
the effect of the on-site U term is to lower the activa- 
tion barrier at z c making the adsorption curve to a fully 
repulsive interaction. Being the chemisorption potential 
energy curve for H on graphene the result of the interplay 
of two an attractive and a repulsive doublet spin states 
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Figure 5: Comparison between GGA (full line) and GGA+U 
(dashed line) potential energy curves along the adsorption 
path. Note how GGA+[7 do not show the correct activation 
barrier arising from the avoided crossing of the two curves 
shown in the right panel of Fig[T] 



their avoided crossing has necessarily to occur. For this 
reason the GGA+f/ results, despite showing the correct 
magnetization, cannot represent the correct H-graphene 



interaction in all its aspects. 



IV. CONCLUSIONS 

The chemisorption of a hydrogen atom onto a flat 
graphene sheet has been studied within semi-local GGA- 
DFT. Contrary to the adiabatic case the substrate sp — 
sp 3 re-hybridization is limited, so adsorbate and sub- 
strate bands become degenerate at a given critical C-H 
distance, z c , where the total magnetization drops to zero. 
This spin transition is due to a fractional spin configu- 
ration: a fictitious effect induced by the self-interaction 
error. To overcome this issue it is possible to use a hy- 
brid functional such as PBEO. The GGA+U approach 
can reproduce the correct magnetization, but leads to 
qualitatively incorrect potential energy curves. 
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